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ABSTRACT 

Using adaptive mesh refinement (AMR) liydrodynamic simulations of vertically stratified hot halo 
gas, we examine the conditions under which clouds can form and condense out of the hot halo medium 
to potentially fuel star formation in the gaseous disk. We find that halo clouds do not develop from 
linear isobaric perturbations. This is a regime where the cooling time is longer than the Brunt- 
Vaisala time, confirming previous linear analysis. We extend the analysis into the nonlinear regime 
by considering mildly or strongly nonlinear perturbations with overdensities up to 100, also varying 
the initial height, the cloud size, and the metallicity of the gas. Here, the result depends on the ratio of 
cooling time to the time required to accelerate the cloud to the sound speed (similar to the dynamical 
time). If the ratio exceeds a critical value near unity, the cloud is accelerated without further cooling 
and gets disrupted by Kelvin- Helmholtz and/or Rayleigh- Taylor instabilities. If it is less than the 
critical value, the cloud cools and condenses before disruption. Accreting gas with overdensities of 
10-20 is expected to be marginally unstable; the cooling fraction will depend on the metallicity, the 
size of the incoming cloud, and the distance to the galaxy. Locally enhanced overdensities within 
cold streams have a higher likelihood of cooling out. Our results have implications on the evolution 
of clouds seeded by cold accretion that are barely resolved in current cosmological hydrodynamic 
simulations and absorption line systems detected in galaxy halos. 

Subject headings: Galaxy: formation — Galaxy: halo — ISM: clouds — methods: numerical 



1. INTRODUCTION 

The Milky Way and other similar-sized star-forming 
disk galaxies will exhaust their cold gas reservoirs within 
a few Gyrs unless it is replenished in some way (e.g, Chi- 
appini et al. 2001; Sommer-Larsen et al. 2003). One 
possible source of this fuel is the direct accretion of cold 
gas from the halo (e.g., Wakker & van Woerden 1997; 
Thilker et al 2004; Oosterloo et al. 2007); however, this 
level of fueling does not appear to be sufficient to balance 
the star formation rate (e.g., Putman 2006; Sancisi et al 
2008; Grcevich k Putman 2009). 

An additional possible source of fuel is gas accretion 
from an extended hot (T ^ lO*' K) halo that is thought 
to surround disk galaxies (e.g., White & Frenk 1991; 
Fukugita & Peebles 2006; Bregman 2007). Ahhough the 
size of this reservoir is not well-constrained, the observa- 
tional indications for its existence are many. These in- 
clude the structure of many high velocity clouds (HVCs; 
Bruens et al. 2000; Peek et al. 2007; Putman et al. 
2011), O VI absorption indicating the HVCs are mov- 
ing through a hot halo medium (Sembach et al. 2003), 
the likelihood that the Magellanic Stream shows signs 
of ram-pressure effects (Weiner & Williams 1996; Put- 
man et al. 2003; Bland-Hawthorn et al. 2007), as well 
as the lack of gas in any dSph galaxy within ~300 kpc 
of the Milky Way (Grcevich & Putman 2009). Observa- 
tional constraints on the temperature and density of the 
hot gas also come from O VII absorption lines in Chan- 
dra and XMM observations of distant quasars and halo 
objects (Wang et al. 2005; Bregman 2007). The exact 
location of this O Vll-absorbing hot gas remains under 
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some debate, but based on the observed line-strength vs. 
Galactic latitude, the majority of the absorption signal is 
construed to lie within ^10 kpc of the plane (Rasmussen 
et al. 2003; Wang et al. 2005; Cohins et al. 2005). On the 
other hand, theoretical models, and some of the obser- 
vational results mentioned above, indicate a very diffuse 
component of the hot gaseous halo extending out be- 
yond 100 kpc (e.g.. Mailer & Bullock 2004; Kaufmann 
et al. 2006; Fukugita & Peebles 2006; Sommer-Larsen 
2006). The potential mass reservoir in this component 
may therefore be substantial (but see Anderson & Breg- 
man 2010). 

Whether or not this hot halo gas can cool and ac- 
crete on to the disk is unclear. One possibility is that 
the gas condenses out directly due to a thermal instabil- 
ity (e.g.. Mailer & Bullock 2004). This idea was exam- 
ined in numerical simulations using the SPH method and 
precisely such condensation was found (Kaufmann et al. 
2006, 2009; Sommer-Larsen 2006) with properties similar 
to those of the observed population of HVCs (Peek et al. 
2008). This is appealing as it explains both the nature 
and source of the HVCs as well as the fuel for ongoing 
star formation in the disk of the Milky Way. 

However, this simple condensation picture has recently 
been challenged. A linear analysis of the thermal insta- 
bility operating in a stable, stratified atmosphere indi- 
cates that the gas is not likely to be unstable unless the 
entropy profile is extremely flat (Binney et al. 2009 and 
references therein). This is unlike the case in an homo- 
geneous medium (e.g. Field 1965; Burkert & Lin 2000). 
The reason is that oscillations driven by the buoyancy 
of the thermally unstable gas will disrupt the gas parcel 
before it can cool. A similar situation has been found to 
exist in clusters of galaxies (e.g. Malagoli et al. 1987; 
Balbus & Soker 1989). If the entropy proflle is flat, or 
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very close to flat, the oscillation timescale can be longer 
than the cooling time and instabilities can occur; how- 
ever, it remains unclear if such profiles actually occur in 
nature, or if they can persist for a long period of time. 

A possible way around this objection is to have the 
perturbations in the gas be nonlinear, such as might be 
found in inflowing filamentary gas which then cools to 
form clouds, as suggested in Keres & Hernquist (2009). 
This is appealing as it links the HVCs and gas replenish- 
ment directly to the relatively cool (T ^ 10^ K), filamen- 
tary fiows found in galaxy formation simulations (Katz 
et al. 2003; Keres et al. 2005; Ocvirk et al. 2008; Brooks 
et al. 2009; Keres et al. 2009). Nonlinear perturbations 
could also result from the leftovers of satellite accretion 
(e.g., Bland-Hawthorn et al. 2007; Grcevich & Putman 
2009). 

In this paper, we numerically investigate the cooling 
of instabilities in a stratified atmosphere. We use pa- 
rameters appropriate for the Milky Way hot halo, and 
investigate both linear and nonlinear instabilities in or- 
der to understand when gas can cool and condense out 
of the hot gas. 

2. METHOD 

We perform simulations with Enzo, an Eulerian hydro- 
dynamics code with adaptive mesh refinement (AMR) 
capability (Bryan & Norman 1997; Norman & Bryan 
1999; Bryan 1999; O'Shea et al. 2004). It solves the 
Euler equations using the piecewise-parabolic method 
(PPM; Colella & Woodward 1984) to handle compress- 
ible flows with shocks. We choose to use PPM rather 
than the Zeus hydro solver, which is also included in 
Enzo, primarily because the former is known to resolve 
instabilities better (e.g., Agertz et al. 2007). 

We set up a stratified isothermal gas at 2.0 x 10^ K 
initially in hydrostatic equilibrium. This value represents 
our best guess for the temperature of the observed hot 
halo gas (e.g., Henley et al. 2010). The computation box 
contains a volume of (1 kpc)^ x (20 kpc), elongated in 
the vertical direction, in order to follow the evolution of 
clouds as they fall towards the midplane (located in the 
middle of our volume). 

The simulations include metallicity-dependent radia- 
tive cooling (Cen et al. 1995) extended down to 10 K 
(Dalgarno & McCray 1972), while evolving a chemical 
gas-phase reaction network of nine species (Abel et al. 
1997). A diffuse form of photoelectric heating (Wolfire 
et al. 1995, 2003) that scales weakly with gas density 
(Joung et al. 2009) is included for gas with T < 10^ K, 
and so is applied only when the cloud cools sufficiently. 

In most of our models, we place up to 3 levels of re- 
finement on top of a base grid of 32^ x 640 cells, re- 
sulting in an effective resolution of 256^ x 5120 cells, 
although we limit refinement only to the top half of the 
box (z > 0). Our refinement criteria are based on gas 
density and cooling time. Specifically, the refinement 
process ensures that the mass contained in a given cell 
is less than 4.9 M© (in the fiducial model) and that the 
cooling time in a cell is always shorter than the sound 
crossing time of that cell. This criterion ensures that 
the cooling time is well resolved. For simplicity, star for- 
mation and energy input from stars are not included in 
our models. We also do not include heat conduction; see 
section 1431 for a discussion of this point. 



The gravitational influence of the disk and dark mat- 
ter halo are included through a static gravitational term, 
which is strictly vertical. The variation of the gravita- 
tional acceleration as a function of height, g{z), is im- 
portant as it determines both the initial density profile 
- assumed to be in hydrostatic equilibrium - and the 
Brunt- Vaisala (B-V) time, i.e. the period of oscillations 
due to buoyancy (defined below). Kuijken & Gilmore 
(1989) derived their estimate of g{z) based on a sam- 
ple of K dwarfs towards the South Galactic pole, which 
is applicable to \z\ < 3 kpc in the solar neighborhood. 
This is employed in our models for \z\ < 3.06 kpc. For 
the layer immediately above it (3.06 kpc < \z\ < 15.2 
kpc), we adopt an expression for the gravitational accel- 
eration in Benjamin & Danly (1997) used to determine 
the distances to falling interstellar clouds in the Galactic 
halo, which in turn was taken from Wolfire et al (1995). 
The value of g{z) within this given range of z is nearly 
constant. Finally, in the uppermost layer, the NFW pro- 
file (Navarro et al. 1997) for a halo with a mass of 2 x 
10^^ Mq and a concentration parameter, c = 12, is used, 
where applicable. Hence, 



a-i- 



a2Z 



g{z) = \ a3tanh(^^^) 



\z\ < 3.06 kpc 



a4fl; \n{l+z/Rs)-z/{Rs+z) 
z'-' ln{l+c)-c/{l+c) 



3.06 kpc < \z\ < 15.2 kpc 
\z\ > 15.2 kpc 

(1) 

Here, ai ^ 1.42 x 10"^ kpc/Myr^, 02 = 5.49 x lO"'^ 
Myr-^ zo = 0.18 kpc, 03 = 9.5 x 10"^ kpc/Myr^, 04 = 
6.97 X 10^9 kpc/Myr^, and Rs = 20 kpc. The three 
curves as well as our choice of g{z) are shown in Figure[TJ 

Figure [2] shows the initial density profile in hydrostatic 
equilibrium with the gravitational field. The overall den- 
sity normalization is not constrained by hydrostatic equi- 
librium; the density is normalized so that n = 10~^ cm"'^ 
at z = 10 kpc (Hsu et al. 2011; Gaensler et al. 2008; 
Rasmussen et al. 2009; Grcevich & Putman 2009; Mc- 
Cammon et al. 2002; Sembach et al. 2003). 

Our model of an initially isothermal, hot gas distribu- 
tion is obviously an over-simplification, since we know 
that near the midplane the mean gas temperature de- 
creases significantly and the mean density is significantly 
higher. However, this simplified model suffices for our 
purpose because we are interested in the early evolution 
of density perturbations, while they are well above the 
disk. Note that including a cooler (e.g., T ~ 10'' K), 
denser midplane layer would have introduced an addi- 
tional difficulty at the interface between the hot and the 
cool layers, as the transition will be smoothed over a few 
cells, and the cooling rate will peak in those interface 
cells, which is unphysical. Ideally, to set up a realistic 
ISM, mechanical (stellar winds and supernova explosions; 
see, e.g., Joung & Mac Low 2006) and radiative feedback 
from star formation should be included, but this is be- 
yond the scope of the present study. 

Periodic boundary conditions are used on the side sur- 
faces perpendicular to the Galactic plane, while outflow 
boundary conditions are used on the upper and lower 
surfaces. We found that the outflow boundary condi- 
tions (more appropriately called 'zero gradient boundary 
conditions') do not lead to a static background medium, 
even if we set up the "active" region of the box (|z| < 10 



Hot Gas Condensation 



3 



kpc in our fiducial model) in exact hydrostatic balance. 
Instead, the background medium gradually developed in- 
flow motions toward the midplane (z = 0), starting from 
the top and bottom boundaries, and after a sufficiently 
long period (> 100 Myr), gas initially located at those 
boundaries moved significantly closer to the midplane. 
The same behavior was seen even after radiative cooling 
was turned off. This is because the outflow boundary 
conditions, which simply copy the variables in the out- 
ermost cell in the active region to the boundary cells, 
neglect the small yet important vertical pressure gradi- 
ent caused by the density drop-off immediately outside 
of the upper and lower boundaries. Note that previous 
models that used outflow boundary conditions in the ver- 
tical direction may have suffered from the same problem, 
although because of the global catastrophic collapse due 
to cooling on shorter timescale, the effect was hidden 
(see, e.g., de Avillez 2000; Joung & Mac Low 2006). As 
a simple remedy for the problem, we set the densities 
in the boundary cells not at the value in the outermost 
cell but at a value slightly lower than it by extending 
the density gradient by one more cell. Doing that was 
enough to stabilize the background medium. 

We introduce a density perturbation with a spherical 
shape (a "cloud") that is overdense with respect to the 
background medium by a prescribed value, which is an 
input parameter of the models. Note that the density in- 
side the cloud is not uniform but is a function of height, 
as the background density decreases in the vertical di- 
rection. The temperature is reduced by the same factor 
so that the cloud is initially in pressure equilibrium with 
the surrounding medium. 

The initial height of the cloud in our fiducial model is 8 
kpc, with a cloud radius of 0.25 kpc, and we use overden- 
sities relative to the mean density oi S — 6p/p — 1.01, 
1.1, 2, 4, 7, 10 and 30 (see Table [1] for parameters for 
all models). In addition, we ran some additional models 
with smaller and larger initial heights, ranging from 2 
kpc to 32 kpc. In runs with our largest heights (16 and 
32 kpc), we increased both the clouds radius and the 
box size to twice and four times the fiducial size in ev- 
ery dimension, respectively, in order to also examine the 
impact of cloud size. The diameter of the cloud is 1/2 of 
the width of the computation box, so the cloud size and 
its mass scale with the box size. The largest height, 32 
kpc, was chosen to be representative of distances where 
cold gaseous streams seed dense clouds aided by hydro- 
dynamic instabilities in the scenario advocated by Keres 
& Hernquist (2009). 

Most simulations were run twice, once with primordial 
abundances (i.e. only H and He cooling), and once with 
a metallicity of 0.3 solar. We ran additional models with 
Z = 0, 0.1 and 0.5 Zq to study the effect of varying 
metallicity. The metallicity of the gas is assumed to be 
uniform throughout the box. 

3. RESULTS 

The simulations we run are described in Table [TJ where 
the naming convention encodes both the height and over- 
density of the cloud. For example, h4-o8 indicates a model 
with a cloud starting at 4 kpc with an overdensity of 8 
relative to the background density at that height. 

3.1. Linear regime 



We begin by examining the results of the model with 
the smallest amplitude in density and temperature (only 
1% of the background values) and compare them with the 
linear predictions made by Binney et al. (2009). They 
predicted that for this background entropy profile, and 
for perturbation amplitudes that are small and approxi- 
mately isotropic, the perturbation should oscillate on the 
B-V time: 



^BV — 27r 



g_ fdT 
T\dz 



7-1 gprup 



7 



-1-1/2 



(2) 



where /i is the mean particle mass and 7 = 5/3 is the ra- 
tio of specific heats. In our fiducial model with an initial 
cloud height of 8 kpc, the cooling time is longer than the 
B-V time by a factor of ^10, and so the cloud (with only 
1% density enhancement) is expected to simply oscillate 
up and down. Half of the time, it is cooler than the 
surrounding gas and half of the time it is hotter, so, to 
first order there should be no net cooling, and the clump 
should not be unstable in the sense of Field (1965). For 
the cooling time, we adopt: 



^cool 



ikT 
2nA{T) 



(3) 



where A(r) is the cooling rate at temperature T. These 
timescales are shown in Figure |3l both as a function of 
clump overdensity and also initial clump height. 

The time sequence in Figure H] shows the z-component 
of the velocity, and confirms that the cloud indeed os- 
cillates vertically around the equilibrium position. The 
grey dotted curve represents a simple sinusoidal wave 
with a period of 527.4 Myr, as predicted by linear anal- 
ysis (eq. [2]). The dashed curve represents the mean ver- 
tical velocity of the cloud in the model without radiative 
cooling, while the solid curve shows the same quantity 
in the model with cooling. The exact shape of the cloud 
is not retained since inside the initially uniform sphere 
various modes having different wavelengths interact with 
one another and distort the waveform. For this reason, 
we plot the mean z-component of the velocity of the os- 
cillation (the area under the vs. x curve, divided by 
the wavelength, here assumed to be the initial diameter 
of the cloucfl) . 

After an initial transient period {t < 100 Myr), the 
mean amplitude of the z-velocity displays a simple sinu- 
soidal shape with a period of 540 Myr, close to the esti- 
mate from the linear analysis, 527.4 Myr. The first half 
period of the cycle, when the cloud is on average below 
the initial position, is slightly shorter than the predicted 
value. On the other hand, the second half period, when it 
is on average above it, is slightly longer by a few percent. 
These differences arise due to the vertical density gradi- 
ent in the background medium as the cloud bobs up and 
down, which results in corresponding changes in the B-V 
times. Regardless of these minor details, this result con- 
firms the prediction by Binney et al. (2009) in the linear 
case (or alternately it can be viewed as a test of the sim- 
ulation method). When the cooling time is significantly 
longer than the B-V time, as in this model (hSol.Ol), 

^ Only the cylindrical volume within Az < 1.5 kpc of the vertical 
line passing through the center of the cloud (Ar < 15.6 pc) was 
used for this analysis. 
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the cloud does not condense or cool but instead simply 
oscillates in the vertical direction due to the interplay of 
gravity and buoyancy. The velocity steadily damps over 
the five oscillation periods shown as the cloud deforms 
and mixes with the background medium. 

3.2. Nonlinear regime 

The results in Binney et al. (2009) were restricted to 
linear analysis and hence only to the onset of the insta- 
bility for very low amplitude perturbations. However, 
recent studies have suggested that nonlinear perturba- 
tions may arise as a natural consequence of the galaxy 
formation process. For example, in the cosmological sim- 
ulation of a Milky Way-like halo by Keres & Hernquist 
(2009), inflowing filamentary gas cools and forms clouds 
at distances of 10-100 kpc from the galaxy. They sug- 
gested this as an important fueling mechanism for contin- 
ued star formation in the disk. It is therefore interesting 
to examine how mildly or highly nonlinear clouds evolve 
over time. 

Figure [5] displays the specific entropy distribution in 
two of our models with initial cloud heights of 8 kpc. 
Each column corresponds to a time slice at (from left to 
right) t = 0, 18, 25, 35, and 65 Myr. With an initial 
overdensity of 4, the cloud does not have sufficient time 
to cool before being disrupted via the Kelvin-Helmholtz 
instability (KHI) (left panels). On the other hand, with 
an initial overdensity of 10 (right panels), the cloud does 
cool and lowers its specific entropy, as can be seen in 
the third column on the right, at i = 25 Myr. After- 
ward, again due to the KHI, the cloud mixes with the 
surrounding hotter gas (last column); however, in this 
paper we focus on the early evolution of the clouds (i.e. 
whether they cool or not) and leave a discussion of their 
ultimate fate to future work. See also Heitsch & Putman 
(2009) for a discussion of the evolution of clouds that are 
fully formed. 

The result shown in Figure [5] implies that there is a 
threshold overdensity, below which the cloud does not 
cool before disruption and above which it does. Lin- 
ear theory does not hold when the cloud overdensity is 
large, as in these cases; however, we can still develop 
an analytic expectation based on time-scale arguments. 
In the analytic case, the criteria for instability can be 
phrased in terms of the ratio of the B-V time to the 
cooling time|3 When tev/icooi < 1, the cloud oscillates 
rather than cooling. In the nonlinear case, the gravita- 
tional force on the cloud dominates the buoyancy force, 
and so we neglect it0 The cloud accelerates relative to 
the background gas and the resulting velocity difference 
leads to surface shear and growth of the KHI. To simplify, 
we assume the KHI results in the rapid destruction of the 
clouds once they reach a significant fraction of the sound 
speed in the background medium, which is ultimately re- 
lated to the depth of the potential well. Therefore, we 
can define an approximate acceleration timescale as 

iacc. = ^=f^)'''. (4) 

g \nmpg^J 

* More precisely, the thermal instability time, although this is 
generally close to the cooling time (see Binney et al. 2009). 

^ Buoyancy will lead to a timescale correction of order , where 
5 is the cloud overdensity. 



Notice that for an isothermal halo, this timescale is di- 
rectly proportional to the B-V time, but is shorter by a 
factor of 7.7 for 7 = 5/3. The criteria for nonlinear cloud 
disruption then becomes approximately icooi/^accoi > 1- 
This simple calculation neglects a number of impor- 
tant effects, such as the time for cloud disruption (which 
we will discuss later), as well as heating due to adiabatic 
compression while the cloud falls. This last effect can 
play a role - for example. Figure |6] shows the time evo- 
lution of the maximum density, minimum temperature, 
and minimum specific entropy in the two models pre- 
sented in Fig. [5l Unlike in the ^ = 10 model (right), 
the (5 = 4 model (left) does not cool significantly before 
disruption at i 60 Myr. The density increase due to 
compression leads to adiabatic heating of the cloud (note 
that the specific entropy is nearly constant in time in the 
left column), which prevents cooling. 

3.2.1. Dependence on cloud height, overdensity, and 
metallicity 

As we have discussed, the key parameter determining 
whether cloud disruption occurs (for nonlinear clouds) is 
the ratio of the cooling time to acceleration time. We 
test this idea using a wide range of simulations, as de- 
scribed in Table [TJ and shown graphically in Figure [71 
The horizontal axis of this figure denotes the overdensity 
of the initially spherical density perturbation (cloud) rel- 
ative to the background medium, while the left and right 
vertical axis show the ratios icooi/^BV and tcooi/^accci, 
respectively. Note that the cloud is set initially in pres- 
sure equilibrium with the surrounding gas and hence the 
temperature inside the cloud is lower by the same over- 
density factor - we use the cloud density and tempera- 
ture to compute the cooling time, while the background 
temperature is used in iev and iaccci- 

Each of the models that we have run are plotted on 
Figure [71 with the symbol indicating the result of the 
simulation. To determine if significant cooling occurs, 
we look at the minimum specific entropy in the cloud, 
and determine the ratio of that minimum entropy to the 
initial entropy. A cross indicates that no significant cool- 
ing occurred (entropy drop less than 20%), a circle is 
used for clouds which clearly cooled (entropy drop larger 
than 50%), while a triangle is used for intermediate cases 
(entropy decrease between 20% and 50%). 

The results define a clear separation between clouds 
that can cool and those that cannot, based on the ratio 
icooi/^accei ~ 1-0 (or cquivalcntly, ^cooi/^BV = 0.12). For 
ratios lower than this critical value, the cloud can cool 
(at least somewhat); while above it, it is disrupted before 
significant cooling can occur. This holds true for a wide 
variety of cloud overdensities and initial heights, and is 
the primary result of this paper. 

We have assumed the gas to have a metallicity of 
0.3Zq, but the metallicity of the hot halo, or inflowing 
clouds, is actually uncertain. High metallicity gas will 
have shorter cooling times and hence be more likely to 
condense out than low metallicity gas, while low metal- 
licity has the opposite effect. To explore this, we have 
repeated all runs with primordial metallicity (shown as 
black symbols in Figure [7]) , and in addition have carried 
out a run with Z — 0.1 Zq and one with Z = 0.5^0 (see 
Table 1 for results) . The change in metallicity affects the 
cooling rate but we see that these models are still con- 



Hot Gas Condensation 



5 



sistent with the relation defined above (i.e. the ratio of 

^cool/^acccl) ■ 

We have also explored the use of the ratio tcooi/^dyn, 
employing tdyn — instead of taccci as the dis- 

ruption timescale. This ratio is quite good in predicting 
whether a given cloud will be able to cool before disrup- 
tion, but does not seem to be as sensitive as tcooi/iaccoi 
for border-line cases. This may not be surprising, as the 
dynamical time is a measure of the time to reach the 
midplane, while the cloud is sensitive to disruption due 
to large velocities. 

3.3. Cloud mass 

The KHI timescale (^kh) is linearly proportional to the 
physical size of the perturbing object (e.g., Nulsen 1982). 
The disruption time is usually found to be a few times 
longer than ^kh, as the instability develops first on the 
smallest scales and must grow to a wavelength compa- 
rable to the size of the object before it can fully disrupt 
the cloud. Therefore, larger clouds might be expected to 
be able to cool before being disrupted. 

To investigate the dependence of cloud disruption time 
on cloud mass, we ran an additional model with hinn = 
32 kpc {h32o30) but with twice the usual cloud radius (2 
kpc instead of 1 kpc). We also double the width of the 
computational volume from 4 x 4 x 40 kpc^ to 8 x 8 x 40 
kpc^ {L8h32o30; see Table [T] for exact parameters). 

In comparing the results of the two runs, which differ 
only in cloud size, we find that the two clouds are accel- 
erated at the same rate, but that the larger cloud takes 
somewhat longer to disrupt. This gives more time for 
cooling to act, and does result in lower entropy and larger 
(by a factor of >1.5) travel distance. This indicates that 
for clouds larger than presented here, tcooi/iaccei may not 
be the only criteria we need to consider when determin- 
ing if a cloud will cool or disrupt. 

The disruption time of the clouds once they cool and 
begin moving through the hot halo is also dependent 
on the mass of the cloud. Heitsch & Putman (2009) 
present a detailed analysis of the disruption time and 
distance traveled for halo clouds similar to the high ve- 
locity clouds at < IC* K that our cooling clouds evolve 
towards. In particular, they found that only clouds with 
masses larger than approximately 10^-^ Mq could make 
it to the midplane from 10 kpc without being disrupted 
(see figure 7 of that paper). In runs where we follow the 
evolution of the cloud after initially cooling, we also find 
the condensations do not travel more than 10 kpc before 
disrupting. This distance is smaller as the mass of the 
cloud decreases - for example, in the hSolO run shown in 
Figure [5] (with a cloud mass of 1.22 x 10^ Mq), disruption 
occurs after about 3-4 kpc. 

3.4. Resolution issues 

Cloud disruption by shear instabilities requires high 
resolution to follow accurately; Klein et al. (1994) 
showed that in the case of a strong shock sweeping over 
a cloud, some aspects of the evolution were only resolved 
with 100 or more cells across the cloud radius, although 
the disruption timescale could be determined with less 
resolution. Our canonical resolution is 64 cells across the 
cloud radius (we denote this as R64). In order to inves- 
tigate the impact of resolution, we have also carried out 



a simulation of the h32ol00 case with R128 (i.e. twice 
the spatial resolution). The results are shown in Figure[8] 
and are very similar for the two different resolution, in- 
dicating that we resolve the cooling and disruption. We 
chose this run for a resolution study because it has the 
largest cell size (since the box is the largest), and hence 
the lowest resolution of any of our standard runs. 

4. DISCUSSION 

What are the implications of these idealized numerical 
models on fueling of the gaseous galactic disks? We dif- 
ferentiate between two sources of cold gas accreting on 
to the disk: (i) gas cooling from the hot halo (T ~ 10^ 
K) via a linear thermal instability, and (ii) gas cooling 
from filamentary flows that have intermediate tempera- 
tures (T 10^ K) and nonlinear overdensities compared 
to the hot halo gas they are flowing through. In sec- 
tion 14.11 we examine the first idea, as it seems to have 
been the motivating mechanism behind the Mailer & Bul- 
lock (2004) method of solving the over-cooling problem 
(although it may not be the only way to achieve multi- 
phase gas). This mechanism is also explicitly examined 
in the SPH simulations of Kaufmann et al. (2006) and 
Sommer-Larsen (2006). The second source (ii) is dis- 
cussed in §4.2 and is relevant to the accretion of cold gas 
from the IGM or satellites. 

4.1. Implications for classical (linear) hot halo cooling 

We find that linear perturbations are generally stable 
in hot galactic halos, in agreement with Binney et al. 
(2009). Therefore, unless the halo or cloud properties 
are quite different from what we model, we do not ex- 
pect cold gas to form via this mechanism. Note that we 
have only examined local cooling instabilities; see Wal- 
ters & Cox (2001) for the possibility of global vertical 
oscillations. 

Before turning to nonlinear instabilities, we examine 
how different conditions could change our results. First, 
as discussed in more detail in Binney et al. (2009), it is, 
in theory, possible to generate linear instabilities via ex- 
tremely oblate perturbations, which increase the pertur- 
bation's inertia and slow the oscillation. This particular 
perturbation configuration seems very unlikely and we 
do not consider it further here. Another possibility is if 
the hot halo has a very flat entropy gradient the buoy- 
ancy is reduced and this puts the ratio icooi/^sv below 
1. Kaufmann et al. (2009) examined this scenario with 
an isolated galaxy formation simulation and initialized 
the gas to have a high uniform entropy, finding that the 
resulting entropy was very flat and that cool clumps did 
form (while in a second simulation which had lower initial 
entropy, they did not). Such high entropy requires a sub- 
stantial amount of non-gravitational heating, more than 
the current generation of cosmological galaxy formation 
simulations produce (e.g. Keres et al. 2009). Finally, 
there are a number of physical effects that we did not 
include in this simulation, including rotation of the halo 
gas, magnetic fields, and thermal conduction. We dis- 
cuss these further in § 14.31 but of all of these, it appears 
that rotation might be the best chance to rescue linear 
perturbations, since some degree of rotation support will 
reduce the buoyancy force and thus the length the B-V 
timescale. However, for this to be important, the rota- 
tional velocity has to be quite large, probably close to 
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Keplerian. 

Since (subject to the above caveats) we do not find 
linear perturbations to condense out of hot galactic gas, 
we turn to the question of why the simulations of Kauf- 
mann et al. (2006) and Sommer-Larsen (2006) did. The 
entropy profiles of these SPH calculations are not given, 
so it is difficult to know if they are significantly flatter 
than we adopt here, but cosmological collapse simula- 
tions generally produce reasonably steep entropy profiles. 
The Sommer-Larsen simulation was cosmological, and so 
did include filamentary flows, and so the perturbations 
might be seeded by these flows (we address this possibil- 
ity in more detail in the next section); however, both sets 
of authors argue that the condensation is coming from 
the hot gas, and Kaufmann et al. (2006) in particular 
suggest that the seeds are relatively small fluctuations 
arising from particle noise. Based on the simulations 
here, we do not expect such fluctuations to cool. One 
other possible reason for this discrepancy is resolution 
- the simulations in this paper have significantly bet- 
ter resolution (several pc) than that in the cosmological 
runs (typically several hundred pc). Another possible 
explanation is the difficulty that SPH codes have resolv- 
ing Kelvin Helmholtz instabiltiies - Agertz et al. (2007) 
found that dense clumps moving through a diffuse back- 
ground were artificially stabilized by particle effects at 
the boundary. This presumably could help explain why 
slightly overdense clouds were not disrupted, and may 
act to suppress the oscillations that stabilize linear over- 
densities. 

4.2. Implications for cloud formation from (nonlinear) 
filamentary flows 

We turn now to the evolution of nonlinear clouds - ini- 
tial perturbations which are well out of the linear regime. 
These might be due to particularly large perturbations 
from dark matter halos, or, more likely, from inflowing 
filamentary gas. In particular, the low-temperature end 
of the Warm Hot Ionized Medium at - 10^ K (WHIM; 
e.g. Dave et al. 2001) will have overdensities of order 
5-30 when flowing into the hot halo. This material is 
sufficiently overdense that it escapes significant heating 
by the virial accretion shock. Clumps of this gas might 
then cool and accrete onto the disk, providing fuel for 
star formation. 

We find that overdensities can cool significantly, pro- 
vided they meet the criteria that the ratio of the cool- 
ing to acceleration time is below about 1 (see Fig. [7]). 
This generally entails substantial overdensities and/or 
low initial cloud heights (since the cooling time drops 
faster than the acceleration time with decreasing height). 
For accretion of WHIM-like gas in filaments at dis- 
tances around 40 kpc, as suggested by Keres & Hern- 
quist (2009), Figures [3] and [7] indicate that overdensities 
of about 10 are required. The overdensities in the fila- 
ment in that paper appear to be similar to this value, 
therefore we conclude that this mechanism for nonlinear 
seeding of cooling instabilities may be a viable way to 
produce cold halo gas. This is one of the key results of 
this paper. We also note that we monitored the material 
cooling and find that the surrounding hot halo medium 
does not cool along with the overdensity. 

The clouds modeled here are initially most similar in 
temperature to the O VI absorbers detected with ultravi- 



olet spectrographs (--10^ K; Thom & Chen 2008; Tripp 
et al. 2008; Prochaska et al. 2011; Narayanan et al. 
2011). The sample of O VI absorbers remains somewhat 
limited, with detections thus far largely beyond 100 kpc 
(in projection) from the nearest galaxy. For most systems 
it remains somewhat unclear if the absorbers are arising 
from collisional or photo-ionization. Recently there has 
been a detection of an O VI absorber at z = 0.35 only 
85 km/s and 95 kpc in projection from a galaxy that is 
found to be a photoionized structure with a length scale 
L « 0.1-1.2 Mpc and log = -4.4 to -4.9 (Thom et al. 
2011). This density is typical of the densities estimated 
for other potentially photoionized O VI absorbers (log 
riH = -4 to < -5; Thom & Chen 2008; Prochaska et al. 
2011). Based on our results, these types of low density 
structures may have difficulty being able to cool. For 
~10^ K flows from the IGM to play a substantial role 
in feeding galaxy halos cold gas, we would expect the 
future absorbers detected by COS to be found closer to 
galaxies and with somewhat higher overdensities. Our re- 
sults may also have implications on the cooling of warmer 
clouds traced by Lya and Mg II absorption line systems, 
but further simulations with initial cloud conditions more 
similar to these should be completed before definitive 
statements can be made. 

For those clouds that do cool before disrupting, they 
usually reach a minimum temperature around 10^ K (al- 
though in some cases T^m < lO'^ K) and continue to 
accelerate. They eventually reach large velocities, and 
are subject to dynamical instabilities which, in our sim- 
ulations, quickly shred the clouds. This can be seen in 
the last panel in Fig. \E\ where the cloud has been dis- 
persed and the averaged specific entropy starts to grow 
with time. This occurs after the cloud has travelled only 
about 10 kpc, and so occurs before it hits the disk. This 
is in agreement with Heitsch & Putman (2009), who also 
find (for similar cloud masses), a maximum distance be- 
fore disruption of about 10 kpc. This seems to imply 
that even if the cloud cools, it will not be able to fuel the 
disk with cold gas, and will instead simply add to the 
hot halo. 

The cloud mass determines whether or not the cloud 
will be disrupted before it reaches the midplane, which 
takes roughly a dynamical time. The larger the mass, 
the longer the disruption length, and so we suspect that 
clouds with even larger masses (than the maximum value 
of 2.07 X 10^ Mq in our models) will travel further, al- 
though this has not specifically been tested. In Keres 
& Hernquist, the typical cloud masses were about 10^ 
Mq , however this may be largely due to their resolution 
limitations. In any case, it would be interesting to fur- 
ther investigate the survival of such large clouds in future 
work as our results indicate the direct accretion of gas 
to the disk from cold flows or stripped satellite gas at 
large radii may be difficult despite the initial ability for 
clumps to cool. 

Finally, we note that there is an interesting regime 
where the cooling time is shorter than the B-V or the 
acceleration time but longer than the dynamical time. 
In this case, the gas will start to cool but will impact 
the disk before it can cool significantly, and so will be 
accreted as warm gas. Such inflowing gas will escape H I 
observations and instead appear as an ionized component 
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with negative velocities (e.g., Haffner et al. 2003). For 
an initial height of 32 kpc (h32) and 0.3 solar metallicity 
gas, this is the case for clouds with initial overdensities 
between ~8 and ^--^15 (see Fig. [3l for primordial gas, the 
range is between ~20 and ^40). Two of our models re- 
side in this regime {h32o30 and L8h32o30 with Z = 0); 
both models indeed show that the clouds stay warm, but 
only in the latter model is it not fully disrupted until 
impacting the midplane. 

4.3. Physical Processes Not Included In the Simulations 

There are a number of physical effects that we do not 
include in our models - we discuss their possible impact 
in this section, starting with magnetic fields, then heat 
conduction, and rotation. 

The impact of the magnetic field on the linear insta- 
bility is subtle - Loewenstein (1990) and Balbus (1991) 
show that, for radial fields, the field can suppress the 
buoyancy oscillations and so enhance the cooling insta- 
bility. However, purely radial fields do not seem a nat- 
ural outcome since they are dynamically weak and even 
mild turbulence will tangle them, leading only to a small 
contribution to the effective isotropic pressure, and neg- 
ligible effect on the instability. 

Magnetic fields may act to stabilize a nonlinear pertur- 
bation against disruption. For the nonlinear case, fields 
may insulate the cloud (e.g. Konz et al. 2002; Dursi 
& Pfrommer 2008), although Stone k Gardiner (2007) 
did not find that a dynamically weak field significantly 
suppressed the growth of Rayleigh- Taylor (R-T) instabil- 
ities. Recently, McClure- Griffiths et al. (2010) measured 
a parallel component of the magnetic field inside an HVC 
associated with the Magellanic System to be greater than 
^6 /iG, which they argue is more than sufficient to sta- 
bilize it against ram pressure stripping. Fields of similar 
strength were not found for 26 other HVCs in this work, 
but this may be partially due to complicated foregrounds. 
It will be an interesting future study to explore similar 
models in a magnetized medium. 

In addition, thermal conduction, which we do not in- 
clude, can act to suppress the cooling of clouds. This is 
treated for linear perturbations in Binney et al (2009), 
but an investigation for nonlinear clouds would also be 
interesting. We expect this would make the clouds even 
more susceptible to disruption, making our results for 
the minimum overdensity for cooling a lower limit. 

Rotation can also impact the stability of linear clouds, 
as centrifugal support will decrease the level of pressure 
support required and therefore lengthen the duration of 
buoyancy oscillations. The level of rotation required for 
this to be important is probably quite large, and so may 
play a role close to the disk, but will be less important 
at large radii. 

Finally, we note that the clouds we model have been 
spherical with a limited range of sizes, and zero initial 
velocity. A more complete investigation of the impact 
of cloud size, initial velocity and cloud shape would be 
useful, but beyond the scope of this paper. 

5. CONCLUSIONS 

Using AMR hydrodynamic simulations of vertically 
stratified hot halo gas, we examine the conditions under 
which clouds can form and condense out to potentially 
fuel star formation in the gaseous disk. 



• We find that halo clouds do not develop from linear 
perturbations in density and temperature. This is 
a regime where the cooling time is longer than the 
B-V time and hence predicted to be stable by Bin- 
ney et al. (2009). Although we did not consider 
this case, if the cooling time were shorter than the 
B-V time at a given height, the entire layer - as op- 
posed to small condensed clouds in the background 
medium - would have cooled and fallen down to 
the midplane, resulting in a global catastrophic col- 
lapse. 

• We extend the theory into the nonlinear regime by 
considering mildly or strongly nonlinear perturba- 
tions with overdensities up to 30 (100 in one model 
with /linit = 32 kpc). Here, the result depends on 
the ratio of cooling time to the time it takes to ac- 
celerate the cloud to a significant fraction of the 
sound speed (i.e. icooi/^accoi)- If this ratio is below 
a critical value, around 1, the cloud can cool signifi- 
cantly before it is disrupted by K-H and/or R-T in- 
stabilities and may reach high enough column den- 
sities to be detectable in HI surveys {Nhi ^ 10^^'^ 
cm^^). 

• We simulate a wide range of cloud overdensities and 
initial heights, and show that the tcooi/^accoi crite- 
ria is a good predictor of whether or not a seeded 
cloud will be able to cool significantly. We show 
this result scales as expected for different metal- 
licities, and is a numerically well-resolved result. 
We also investigate cloud size and find that, once 
formed, large clouds take longer to disrupt. 

• Nonlinear perturbations, such as those seeded by 
filamentary flows from cosmological simulations 
(e.g. Keres and Hernquist 2009), appear to be 
right at the critical threshold and therefore proba- 
bly can cool before being disrupted. However, we 
find (in agreement with Heitsch & Putman (2009)) 
that clouds with M < 10^ M0 are disrupted within 
10 kpc, and suspect that only larger mass clouds 
can potentially survive the trip to the disk. 
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Figure 1. Vertical gravitational acceleration adopted in our models. The thick light grey curve represents our choice of g(z) at different 
ranges of heights when each analytic expression is thought to be valid; see equation In the figure, NFW indicates the (radial) acceleration 
from Navarro et al (1997), KG is the vertical acceleration from Kuijken & Gilmore (1989), appropriate for z < 3 kpc, and the line marked 
BD is from Benjamin & Danly (1997). 
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Figure 2. Initial density profiles in hydrostatic equilibrium with the gravitational field in equation jlj. The density is normalized to 10~* 
cm^'^ at 2 = 10 kpc. The vertical lines at 2 = 3.06 and 15.2 kpc indicate where the analytic expression for g{z) changes. The solid curve is 
the density profile for T^aio = 2 X 10® K and fi = 0.61, as adopted in our models. The dashed curve shows how the profile would change if 
we adopted a lower (but still consistent with observational constraints) value for Thaio of 1 X 10® K. For the dotted curve, Thaio = 2 X 10® 
K is again assumed but /i = 1.2. 
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Figure 3. The solid lines show cooling times as a function of the initial overdensity for a range of initial cloud heights from 1 to 32 kpc, 
from bottom to top. For these curves, we assume a metallicity of 0.3 solar; the thick dotted line in dark grey shows the cooling times 
for 0.1 solar metallicity gas at 32 kpc, while the other thick dotted line in light grey shows those for a background medium that is only 
10% as dense as in the fiducial model. Also shown, as dashed lines, are the B-V (left) and dynamical times (right) for the same region 
of initial cloud heights. The exact shape of the cooling times are dictated by the cooling function, A(T), as the initial temperature is 
inversely proportional to the overdensity factor (note that pressure equilibrium was assumed). Hence the cooling times have their minima 
at Pci/Pbg ~ 200, which corresponds to an initial cloud temperature of T^aio / 200 ~ 10^ K. The B-V and dynamical times do not depend 
on the cloud overdensity as they are determined by properties of the background medium and the gravitational acceleration at a given 
height. The cloud will cool if tcoo; < ^BVy otherwise it will oscillate vertically. 
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Figure 4. The vertical velocity of the cloud as a function of time in model hSol.Ol with 5 = 0.01. The grey dotted curve represents a 
perfect sinusoidal save with a period of 527.4 Myr (represented by the vertical dot-dashed line) as predicted by linear analysis. The dark 
dashed curve shows the mean vertical velocity of the cloud in the model without radiative cooling. The dark solid curve shows the same 
quantity in the model with cooling. In the model with cooling, the cloud slowly sinks due to cooling of the atmosphere below the cloud. 
To make the plot clearer, we have removed this nearly constant negative z-velocity from the plotted curves. It was measured by averaging 
the vertical velocity well above and below the cloud, and is shown as the light grey solid and dashed lines for the cooling and no-cooling 
cases, respectively. 
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LOG Specific Entropy {< MpcVwsun^''^) 
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Figure 5. Specific entropy distribution in two of our models with the initial cloud height of 8 kpc. Each column corresponds to a time 
slice at (from left to right) t = 0, 18, 25, 35, and 65 Myr. (a) With an initial overdensity of 4, the cloud does not have sufficient time to 
cool before being disrupted via KHI. (b) With an initial overdensity of 10, the cloud cools and condenses, lowering the specific entropy 
of the cloud (see the third column). Afterward, again due to KHI, the cloud mixes with the surrounding hotter gas (last column). The 
colorbar shown at the bottom applies to both models. 
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Figure 6. Time evolution of maximum density (dotted), minimum temperature (dashed), and minimum specific entropy (solid) in the two 
models presented in Fig. (5] showing the fractional change from their respective initial values. (Note that the vertical ranges are different 
in the two plots.) Unlike in the 5 = 10 model (right), the 5 = 4 model (left) does not cool significantly before disruption at t f» 60 Myr. 
The density increase due to compression leads to adiabatic heating of the cloud (note that the specific entropy is nearly constant in time in 
the left column), which slows cooling. For these figures, cells with the minimum 3% in terms of specific entropy were used, with averages 
weighted by mass. 
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Figure 7. We plot the models we have run in the space of overdensity and the ratio of the cooling time to the B-V time (left) scale, or 
to taccel ('''ight) scale. Horizontal dashed lines show unity ratios for both scales. Crosses indicates runs that clearly did not cool (entropy 
drop by less than 20%), triangles are marginal cases (entropy drop between 20% and 50%, while circles are runs which show clear cooling 
(central entropy drop larger than 50%). Black symbols assume primordial abundances, while green symbols are used for metallicity of 
Z = O.SZq. 




Figure 8. A comparison of the maximum density and minimum temperature for the h32ol00 run at two different resolutions: 64 cells 
across the cloud radius (R64, black) and 128 cells across the cloud radius (R128, red). Noticeable differences are seen only after i 60 
Myr, when the cloud is gradually disrupted due to the KHI, which is not the main focus of this paper. 
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Table 1 

Parameters for the simulations 



model I /linit (kpc) "bg (cm ^) 5 Z {Zq) Rcioud (kpc) ^cioud (M©) Cooled? 



ft.2ol.01 


2 


2.35 X 10~* 


1.01 


0, 0.3 


0.25 


2.34 X 10^ 


N, N 


ft2ol.l 






1.1 






2.54 X 10^ 


N, N 


h2o2 






2 






4.62 X 10^ 


N, N 


h2oi 






4 






9.25 X 10^ 


N P 


h2o7 






7 






1.62 X 10^ 


P P 


h2oW 






10 






2.31 X 10^ 


Y Y 


/i2o30 






30 






6.94 X 10^ 


Y Y 


Mol.Ol 


4 


1.91 X 10~* 


1.01 


0, 0.3 


0.25 


1.90 X 10^ 


N, N 


Mol.l 






1.1 






2.07 X 10^ 


N, N 


Mo2 






2 






3.76 X 10^ 


N, N 


hAoi 






4 






7.52 X 10^ 


N, N 


hAo7 






7 






1.32 X 10^ 


N, N 


hioW 






10 






1.88 X 10^ 


Y Y 


h4o30 






30 






5 64 X 10^ 


Y Y 


hf^n^ 01 


8 


1.24 X 10~* 


1.01 


0, 0.3 


0.25 


j-.^tj A, -Lyj 


N N 


hSol.l 






1.1 






1.34 X 10"^ 


N N 








9 








IN , In 


hSo4 






4 






4.88 X 10^ 


N, N 


h8o7 






7 






8.54 X 10^ 


N, Y 


hSolO 






10 






1.22 X 103 


P, Y 


h8o30 






30 






3.66 X 10^ 


Y, Y 


fcl6o30 


16 


5.20 X lO"-' 


30 




0.5 


1.23 X 10* 


P, Y 


/i32ol0 


32 


1.37 X 10-5 


10 




1 


8.63 X 10'^ 


N, Y 


fe32o30 






30 






2.59 X 10"' 


P, Y 


fe32ol00 






100 






8.63 X 10"* 


Y, Y 


?i32o30z0.1 






30 


0.1 




2.59 X 10-* 


Y 


/i32o30z0.5 








0.5 






Y 


h32o30rh 








0.3 


0.5 


3.24 X 10-'' 


Y 


h32o30rd 










2 


2.07 X 10^ 


Y 


L8/i32o30 








0, 0.3 






P, Y 



Note. — The columns show the model name, initial cloud height, number density of the background medium at /linit, cloud overdensity, 
assumed metallicity of gas, initial cloud radius, and whether the cloud cooled significantly (Y), partly (P) or not (N). Blank indicates the 
value remains the same as that in the above row. 



